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Abstract 

We present the results of searches for dipolar-type anisotropics in different energy ranges 
above 2.5 x 10 17 eV with the surface detector array of the Pierre Auger Observatory, report- 
ing on both the phase and the amplitude measurements of the first harmonic modulation 
in the right-ascension distribution. Upper limits on the amplitudes are obtained, which 
provide the most stringent bounds at present, being below 2% at 99% C.L. for EeV ener- 
gies. We also compare our results to those of previous experiments as well as with some 
theoretical expectations. 



1. Introduction 

The large-scale distribution of the arrival directions of Ultra-High Energy Cosmic Rays 
(UHECRs) is, together with the spectrum and the mass composition, an important observ- 
able in attempts to understand their nature and origin. The ankle, a hardening of the energy 
spectrum of UHECRs located at E 4 EeV , where 1 EeV = 10 18 eV, is presumed 
to be either the signature of the transition from galactic to extragalactic UHECRs or 
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the distortion of a proton-dominated extragalactic spectrum due to e ± pair production of 
protons with the photons of the Cosmic Microwave Background (CMB) [a [ZD. If cosmic 
rays with energies below the ankle have a galactic origin, their escape from the Galaxy 
might generate a dipolar large-scale pattern as seen from the Earth. The amplitude of such 
a pattern is difficult to predict, as it depends on the assumed galactic magnetic field and the 
charges of the particles as well as the distribution of sources. Some estimates, in which the 
galactic cosmic rays are mostly heavy, show that anisotropics at the level of a few percent 
are nevertheless expected in the EeV range Jsl Even for isotropic extragalactic cosmic 
rays, a dipole anisotropy may exist due to our motion with respect to the frame of extra- 
galactic isotropy. This Compton-Getting effect ft 1 QTI has been measured with cosmic rays of 
much lower energy at the solar frequency 111 1 L Il2ll as a result of our motion relative to the 
frame in which they have no bulk motion. 

Since January 2004, the surface detector (SD) array of the Pierre Auger Observatory 
has collected a large amount of data. The statistics accumulated in the 1 EeV energy range 
allows one to be sensitive to intrinsic anisotropics with amplitudes down to the 1% level. 
This requires determination of the exposure of the sky at a corresponding accuracy (see 
Section [3]) as well as control of the systematic uncertainty of the variations in the count- 
ing rate of events induced by the changes of the atmospheric conditions (see Section HI). 
After carefully correcting these experimental effects, we present in Section |5] searches for 
first harmonic modulations in right-ascension based on the classical Rayleigh analysis |fl3ll 
slightly modified to account for the small variations of the exposure with right ascension. 

Below E 1 EeV, the detection efficiency of the array depends on zenith angle and 
composition, which amplifies detector-dependent variations in the counting rate. Conse- 
quently, our results below 1 EeV are derived using simple event counting rate differences 
between Eastward and Westward directions Q14[|. That technique using relative rates allows 
a search for anisotropy in right ascension without requiring any evaluation of the detection 
efficiency. 

From the results presented in this work, we derive in Section |6] upper limits on modu- 
lations in right-ascension of UHECRs and discuss some of their implications. 



2. The Pierre Auger Observatory and the data set 

The southern site of the Pierre Auger Observatory [15J is located in Malargiie, Ar- 
gentina, at latitude 35.2° S, longitude 69.5° W and mean altitude 1400 meters above sea 
level. Two complementary techniques are used to detect extensive air showers initiated by 
UHECRs : a surface detector array and & fluorescence detector. The SD array consists of 
1660 water-Cherenkov detectors covering an area of about 3000 km 2 on a triangular grid 
with 1.5 km spacing, allowing electrons, photons and muons in air showers to be sam- 
pled at ground level with a duty cycle of almost 100%. In addition, the atmosphere above 
the SD array is observed during clear, dark nights by 24 optical telescopes grouped in 4 
buildings. These detectors observe the longitudinal profile of air showers by detecting the 
fluorescence light emitted by nitrogen molecules excited by the cascade. 

The data set analysed here consists of events recorded by the surface detector from 1 
January 2004 to 3 1 December 2009. During this time, the size of the Observatory increased 



6 



from 154 to 1660 surface detector stations. We consider in the present analysis eventsQ with 
reconstructed zenith angles smaller than 60° and satisfying a fiducial cut requiring that the 
six neighbouring detectors in the hexagon surrounding the detector with the highest signal 
were active when the event was recorded. Throughout this article, based on this fiducial 
cut, any active detector with six active neighbours will be defined as an unitary cell \\6\ . It 
ensures both a good quality of event reconstruction and a robust estimation of the exposure 
of the SD array, which is then obtained in a purely geometrical way. The analysis reported 
here is restricted to selected periods to eliminate unavoidable problems associated to the 
construction phase, typically in the data acquisition and the communication system or due 
to hardware instabilities [16]. These cuts restrict the duty cycle to - 85%. Above the 
energy at which the detection efficiency saturates, 3 EeV [16], the exposure of the SD array 
is 16,323 km 2 sr yr for six years used in this analysis. 

The event direction is determined from a fit to the arrival times of the shower front at 
the SD. The precision achieved in this reconstruction depends upon the accuracy on the 
GPS clock resolution and on the fluctuations in the time of arrival of the first particle [11711 . 
The angular resolution is defined as the angular aperture around the arrival directions of 
cosmic rays within which 68% of the showers are reconstructed. At the lowest observed 
energies, events trigger as few as three surface detectors. The angular resolution of events 
having such a low multiplicity is contained within 2.2°, which is quite sufficient to perform 
searches for large-scale patterns in arrival directions, and reaches ~ 1° for events with 
multiplicities larger than five [|18||. 

The energy of each event is determined in a two-step procedure. First, using the con- 
stant intensity cut method, the shower size at a reference distance of 1000 m, 5(1000), is 
converted to the value S 38 ° that would have been expected had the shower arrived at a zenith 
angle 38°. Then, S 38 ° is converted to energy using a calibration curve based on the fluores- 
cence telescope measurements Ill9ll . The uncertainty in S 38 ° resulting from the adjustment 
of the shower size, the conversion to a reference angle, the fluctuations from shower-to- 
shower and the calibration curve amounts to about 15%. The absolute energy scale is given 
by the fluorescence measurements and has a systematic uncertainty of 22% II 1911 . 



3. The exposure of the surface detector 

The instantaneous exposure co(t, 6, (p, .S 38 =) of the SD array at the time t as a function of 
the incident zenith and azimuth^ angles (6, <p) and shower size S 38 = is given by : 

(i)(t,6,<p,S3&°) =n ce u(OX«cellCOS0Xe(S3 8 o,0,0), (1) 

where a ce n cos 6 is the projected surface of a unitary cell under the incidence zenith angle 6, 
"ceii(0 is the number of unitary cells at time t, and e(S 38 °, 6, (p) is the directional detection 
efficiency at size parameter S 38 ° under incidence angles (0, (p). The conversion from S 38 = to 
the energy E, which accounts for the changes of atmospheric conditions, will be presented 
in the next section. 



'A comprehensive description of the identification of shower candidates detected at the SD array of the 
Pierre Auger Observatory is given in reference [11 61 . 

2 The angle (f> is the azimuth relative to the East direction, measured counterclockwise. 
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Figure 1: Relative variation of the integrated number of unitary cells as a function of the solar hour 
of the day in UTC (left panel), and as a function of the local sidereal time (right panel). 



The number of unitary cells n ce \\(t) is recorded every second using the trigger system of 
the Observatory and reflects the array growth as well as the dead periods of each detector. It 
ranges from ^ 60 (at the begining of the data taking in 2004) to ^ 1200 (from the middle of 
2008). From Eqn. \T\ it is apparent that rc cell (f) is the only time-dependent quantity entering 
in the definition of the instantaneous exposure, modulating within any integrated solid angle 
the expected number of events as a function of time. For any periodicity T, the total number 
of unitary cells N ce n(t) as a function of time t within a period and summed over all periods, 
and its associated relative variations AN ce n(t) are obtained from : 



N cen (t) 
(N cen (t)Y 



(2) 



with<AU(0> = 1/r £dtN csll (t). 

A genuine dipolar anisotropy in the right ascension distribution of the events induces 
a modulation in the distribution of the time of arrival of events with a period equal to one 
sidereal day. A sidereal day indeed corresponds to the time it takes for the Earth to complete 
one rotation relative to the vernal equinox. It is approximately T sid = 23 hours, 56 minutes, 
4.091 seconds. Throughout this article, we denote by a the local sidereal time and express 
it in hours or in radians, as appropriate. For practical reasons, a is chosen so that it is 
always equal to the right ascension of the zenith at the center of the array. 

On the other hand, a dipolar modulation of experimental origin in the distribution of 
the time of arrival of events with a period equal to one solar day may induce a spurious 
dipolar anisotropy in the right ascension distribution of the events. Hence, it is essential 
to control AN ce u(t) to account for the variation of the exposure in different directions. We 
show AA^ ce ii(0 m Fig. Q] in 360 bins of 4 minutes at these two time scales of particular 
interest : the solar one T = T so i = 24 hs (left panel), and the sidereal one T = T sie j (right 
panel). A clear diurnal variation is apparent on the solar time scale showing an almost 
dipolar modulation with an amplitude of ^ 2.5%. This is due to both the working times 
of the construction phase of the detector and to the outage of some batteries of the surface 
detector stations during nights. When averaged over 6 full years, this modulation is almost 
totally smoothed out on the sidereal time scale as seen in the right panel of Fig. \T\ This 
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distribution will be used in section l5T.ll to weight the events when estimating the Rayleigh 
parameters. 

From the instantaneous exposure, it is straightforward to compute the integrated expo- 
sure either in local coordinates co(6, (p, a ) by replacing n ce n(i) by AA^ ce ii(a°) in Eqn. [H or 
in celestial coordinates co(a, 6) by expressing the zenith angle 8 in terms of the equatorial 
right ascension a and declination 5 through : 

cos 6 = sin ^ site sin 6 + cos 4ite cos 6 cos (a - a ), (3) 

(where i site is the Earth latitude of the site) and then by integrating Eqn. [Hover time. Be- 
sides, let us also mention that to account for the spatial extension of the surface detector 
array making the latitude of the site 4ite varying by ^ 0.4°, the celestial coordinates (a, 6) 
of the events are calculated by transporting the showers to the "center" of the Observatory 
site. 



4. Influence of the weather effects 

Changes in the atmospheric pressure P and air density p have been shown to affect 
the development of extensive air showers detected by the surface detector array and these 
changes are reflected in the temporal variations of shower size at a fixed energy [20]. To 
eliminate these variations, the procedure used to convert the observed signal into energy 
needs to account for these atmospheric effects. This is performed by relating the signal 
at 1 km from the core, 5(1000), measured at the actual density p and pressure P, to the 
one 5(1000) that would have been measured at reference values po and Po, chosen as the 
average values at Malargiie, i.e. p = 1.06 kg m~ 3 and P = 862 hPa jioll : 

5 (1000) = [l - a P (6)(P - Po) - a p (9)(p d - p„) - /3 p (9)(p - p d )] 5 (1000), (4) 

where p c [ is the average daily density at the time the event was recorded. The mea- 
sured coefficients a p = (-0.80 + 0.02) kg'm 3 , fi p = (-0.21 + 0.02) kg-W and a P = 
(-1.1 +0.1) 10" 3 hPa 1 reflect respectively the impact of the variation of air density (and 
thus temperature) at long and short time scales, and of the variation of pressure on the 
shower sizes [20]. It is worth pointing out that air density coefficients are here predominant 
relative to the pressure one. The zenithal dependences of these parameters, that we use in 
the following, were also studied in reference [20]. It is this reference signal 5(1000) which 
has to be converted, using the constant intensity cut method, to the signal size 5 38° an d 
finally to energy. For convenience, we denote hereafter the uncorrected (corrected) shower 
size 5 3 8 = (5 38°) simply by 5 (5). 

Carrying out such energy corrections is important for the study of large scale anisotro- 
pics. Above 3 EeV, the rate of events R per unit time above a given uncorrected size 
threshold 5,/,, and in a given zenith angle bin, is modulated by changes of atmospheric 
conditions : 



R(>S th )oc dS-^r — oc l+( 7s -l)^A£ dSS 
Js th dS ">j L J J s 



~ ys , (5) 

'Sth 
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Figure 2: Relative modulation of the rates above a given corrected signal size due to the variations of 
the detection efficiency under changes of atmospheric conditions, relatively to the factor modulating 
the rate of events above the corresponding uncorrected signal size in units of S 1/2- 



where hereafter £ genetically denotes P, p or p^, and where we have adopted for the differ- 
ential flux dJ/dS a power law with spectral index y s . Hence, under changes of atmospheric 
parameters A£, the following relative change in the rate of events is expected : 

l dR(>S th ) 

R-^f- * to" 1 )"*- < 6) 

Over a whole year, this spurious modulation is partially compensated in sidereal time, 
though not in solar time. In addition, a seasonal variation of the modulation of the daily 
counting rate induces sidebands at both the sidereal and anti-sidereal frequencies. This 
may lead to misleading measures of anisotropy if the amplitude of the sidebands signifi- 
cantly stands out above the background noise [21]. Correcting energies for weather effects, 
the net correction in the first harmonic amplitude in sidereal time turns out to be only of 
- 0.2% for energy thresholds greater than 3 EeV, thanks to large cancellations taking place 
when considering the large time period used in this study. 

In addition to the energy determination, weather effects can also affect the detection 
efficiencies for showers with energies below 3 EeV, for the detection of which the surface 
array is not fully efficient. Changes of the shower signal size due to changes of weather 
conditions A£ imply that showers are detected with the efficiency associated to the observed 
signal size S , which is related at first order to the one associated to the corrected signal size 
through : 



e(S)^e(S) + (S-S) d(iS} 



dS 



s=s 



dS 



(7) 

5=5 



3 The anti-sidereal time is a fictitious time scale symmetrical to the sidereal one with respect to the solar 
one and that reflects seasonal influences [21]. It corresponds to a fictitious year of ^ 364 days. 
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where we have made use of Eqn. HI The second term modulates the observed rate of 
events, even after the correction of the signal sizes. Indeed, the rate of events R above a 
given corrected signal size threshold S t h is now the integration of the cosmic ray spectrum 
weighted by the corresponding detection efficiency expressed in terms of the observed 
signal size S 

de(S) 



/?(> S,„) « f dS 

J Sth 



e(S) + a^S 



dS 



5=5 



dJ 
dS' 



(8) 



Hence, the relative change in the rate of events under changes in the atmosphere becomes : 



1 dR(> S th ) 
R Ik 



R Js th dS 



de(S) 



dS 



(9) 



5=5 



which, after an integration by parts and at first order in a f , leads to : 



1 dR(> S th ) 
R d£ 



(y s - 1)0£ 



1 - 



e(S th ) CdS§ 

J J ih 



-ys 



I d§e(S)S 



■ys 



(10) 



The expression in brackets gives the additional modulations (in units of the weather effect 
modulation (y s - \)a^ when the detection efficiency is saturated) due to the variation of 
the detection efficiency. Note that this expression is less than 1 for any rising function e 
satisfying < e(S) < 1, and reduces to 0, as expected, when e(S t h) = 1- As a typical 
example, we show in Fig. [2] the expected modulation amplitude as a function of S by 
adopting a reasonable detection efficiency function of the form e(S) = S 3 /[S 3 + 5^ /2 ] 
where the value of S i/ 2 is such that e(S i/ 2 ) = 0.5. This relative amplitude is about 0.3 for 
S = S 1/2, showing that for this signal size threshold the remaining modulation of the rate 
of events after the signal size corrections is about 0.3 x (y s - \)a^. The value of S 1/2 being 
such that it corresponds to - 0.7 EeV in terms of energy, it turns out that within the current 
statistics the Rayleigh analysis of arrival directions can be performed down to threshold 
energies of 1 EeV by only correcting the energy assignments. 



5. Analysis methods and results 

5.1. Overview of the analyses 

The distribution in right ascension of the flux of CRs arriving at a detector can be 
characterised by the amplitudes and phases of its Fourier expansion, 1(a) = I (l+r cos(ar - 
(p) + r' cos(2(a - (£')) + ...). Our aim is to determine the first harmonic amplitude r and 
its phase (p. To account for the non-uniform exposure of the SD array, we perform two 
different analyses. 

5.1.1. Rayleigh analysis weighted by exposure 

Above 1 EeV, we search for the first harmonic modulation in right ascension by apply- 



ing the classical Rayleigh formalism [13] slightly modified to account for the non-uniform 



exposure to different parts of the sky. This is achieved by weighting each event with a 
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factor inversely proportional to the integrated number of unitary cells at the local sidereal 



time of the event (given by the right panel histogram of Fig. [TJ) 11221 12311 



2 N 2 N 

a = — 2_j w i cos a u b = — 2^ Wi sin a u (11) 

i=\ (=1 

where the sum runs over the number of events N in the considered energy range, the weights 
are given by w, = [A/V ce ii(a^)] -i and the normalisation factor is N = Yli=\ W;. The estimated 
amplitude r and phase <p are then given by : 

r = Va 2 + b 2 , ip = arctan -. (12) 

a 

As the deviations from an uniform right ascension exposure are small, the probability P(> 
r) that an amplitude equal or larger than r arises from an isotropic distribution can be 
approximated by the cumulative distribution function of the Rayleigh distribution P(> r) = 
exp (-&o), where k = Nr 2 /4. 

5.1.2. East-West method 

Below 1 EeV, due to the variations of the event counting rate arising from Eqn. [lOl 



we adopt the differential East-West method Ill4ll . Since the instantaneous exposure of the 



detector for Eastward and Westward events is the same0, with both sectors being equally 
affected by the instabilities of the detector and the weather conditions, the difference be- 
tween the event counting rate measured from the East sector, Ie(ch ), and the West sector, 
I w (a°), allows us to remove at first order the direction independent effects of experimental 
origin without applying any correction, though at the cost of a reduced sensitivity. Mean- 
while, this counting difference is directly related to the right ascension modulation r by 
(see Appendix) : 

n n N 2 (sin 9) „ 
I E (a°) - I w (a°) = r sin(a° - <p). (13) 

2n 7T(COS 0) 

The amplitude r and phase (p can thus be calculated from the arrival times of each set of Af 
events using the standard first harmonic analysis N 1 3M slightly modified to account for the 
subtraction of the Western sector to the Eastern one. The Fourier coefficients a EW and b EW 
are thus defined by : 

2 n 2 N 

a EW = JjYj C ° S ^ + ^' bEW = A? ^ ^ + ( 14 ) 

(' = 1 (' = 1 

where ^ equals if the event is coming from the East or n if coming from the West (so 
as to effectively subtract the events from the West direction). This allows us to recover the 
amplitude r and the phase (f EW from 



ew + b E\v> 'Pew = arctan | - — |. (15) 



b 

®EW 



The global tilt of the array of ^ 0.2°, that makes it slightly asymmetric, is here negligible - see sub- 
section 



12 



4 



CD 
"D 
=i 

"Q. 

E 
< 



No correction 
Energy Correction 
' + Exposure correction 




363.5 364 364.5 365 365.5 366 366.5 367 
Frequency [cycles/year] 



Figure 3: Amplitude of the Fourier modes as a function of the frequency above 1 EeV. Thin dotted 
curve : before correction of energies and exposure. Dashed curve : after correction of energies but 
before correction of exposure. Thick curve : After correction of energies and exposure. Dashed 
vertical lines from left to right : anti-sidereal, solar and sidereal frequencies. 



Note however that ip^w, being the phase corresponding to the maximum in the differential 
of the East and West fluxes, is related to <p in Eqn. [T2~l through <p = (p EW + n/2. As in the 
previous analysis, the probability P(> r) that an amplitude equal or larger than r arises 
from an isotropic distribution is obtained by the cumulative distribution function of the 
Rayleigh distribution P{> r) = exp(-fc ( f w ), where k™ = (2 (sinfl) /tt(cos 5)) 2 x Nr 2 /4. 
For the values of (sin 6) and (cos 8) of the events used in this analysis, the first factor in the 
expression for k^ w is 0.22. Then, comparing it with the expression for k in the standard 
Rayleigh analysis, it is seen that approximately four times more events are needed in the 
East- West method to attain the same sensitivity to a given amplitude r. 

5.2. Analysis of solar, anti-sidereal, and random frequencies 

The amplitude r corresponds to the value of the Fourier transform of the arrival time 
distribution of the events at the sidereal frequency. This can be generalised to other fre- 
quencies by performing the Fourier transform of the modified time distribution Il24ll : 

a^^ + ai-a*. (16) 

1 sid 

Such a generalisation is helpful for examining an eventual residual spurious modulation 
after applying the Rayleigh analysis after the corrections discussed in Sections [3] and |4] 
The amplitude of the Fourier modes when considering all events above 1 EeV are shown in 
Fig.[3]as a function of the frequency in a window centered on the solar one (indicated by the 
dashed line at 365.25 cycles/year). The thin dotted curve is obtained without accounting 
for the variations of the exposure and without accounting for the weather effects. The large 
period of time analysed here, over 6 years, allows us to resolve the frequencies at the level 
of ^ 1/6 cycles/year. This induces a large decoupling of the frequencies separated by more 
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Figure 4: Left : Relative variation of the integrated number of unitary cells as a function of the time 
fioo, where the time scale is such that corresponding frequency is 100 cycles/year. Right : Rayleigh 
analysis above 1 EeV for 1600 random frequencies ranging from 100 to 500 cycles/year. Thin 
histogram : analysis without accounting for the exposure variations. Thick histogram : analysis 
accounting for the exposure variations. Smooth curve : expected Rayleigh distribution. 



than this resolution 112411 . In particular, as the resolution is less than the difference between 
the solar and the (anti-)sidereal frequencies (which is of 1 cycle/year), this explains why 
the large spurious modulations standing out from the background noise around the solar 
frequency are largely averaged out at both the sidereal and anti- sidereal frequencies even 
without applying any correction. The impact of the correction of the energies discussed 
in Section |4] is evidenced by the dashed curve, which shows a reduction of - 30% of the 
spurious modulations within the resolved solar peak. In addition, when accounting also for 
the exposure variation at each frequency, the solar peak is reduced at a level close to the 
statistical noise, as evidenced by the thick curve. Results at the solar and the anti-sidereal 
frequencies are collected in Tab. [D 





7" so lar[%] 


P(> r S0 lar)[%] 




P(> r a nti-sid)[%] 


no correction 


3.7 


* 2 10" 37 


0.36 


43 


energy corrections 


2.9 


* 4 10~ 23 


0.15 


85 


+exposure correction 


0.96 


0.2 


0.49 


19 



Table 1: Amplitude and corresponding probability to get a larger amplitude from an isotropic distri- 
bution at both the solar and the anti-sidereal frequencies for events with energies > 1 EeV. 



To provide further evidence of the relevance of the corrections introduced to account for 
the non-uniform exposure, it is worth analysing on a statistical basis the behaviour of the 
reconstructed amplitudes at different frequencies (besides the anti- sidereal/solar/sidereal 
ones). In particular, as the number of unitary cells rc ce n has increased from 60 to 1200 
over the 6 years of data taking, an automatic increase of the variations of A/V cell (?) is ex- 
pected at large time periods. This expectation is illustrated in the left panel of Fig.ffl which 
is similar to Fig. [T]but at a time periodicity T - 87.5 h, corresponding to a low frequency 
of 100 cycles/year. The size of the modulation is of the order of the one observed in Fig. Q] 
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Figure 5: Top : Amplitude of the first harmonic as a function of energy. The dashed line indicates the 
99% C.L. upper bound on the amplitudes that could result from fluctuations of an isotropic distri- 
bution. Bottom : Corresponding probabilities to get at least the same amplitude from an underlying 
isotropic distribution. 



at the solar frequency. In the right panel of Fig. |4] the results of the Rayleigh analysis 
applied above 1 EeV to 1600 random frequencies ranging from 100 to 500 cycles/year are 
shown by histograming the reconstructed amplitudes. The thin one is obtained without ac- 
counting for the variations of the exposure : it clearly deviates from the expected Rayleigh 
distribution displayed in the same graph. Once the exposure variations are accounted for 
through the weighting procedure, the thick histogram is obtained, now in agreement with 
the expected distribution. Note that in both cases the energies are corrected for weather 
effects, but the impact of these effects is marginal when considering such random frequen- 
cies. This provides additional support that the variations of the counting rate induced by 
the variations of the exposure are under control through the monitoring of AN ce \\(t). 
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Figure 6: Phase of the first harmonic as a function of energy. The dashed line, resulting from an 
empirical fit, is used in the likelihood ratio test (see text). 



5.3. Results at the sidereal frequency in independent energy bins 

To perform first harmonic analyses as a function of energy, the choice of the size of 
the energy bins, although arbitrary, is important to avoid the dilution of a genuine signal 
with the background noise. In addition, the inclusion of intervals whose width is below 
the energy resolution or with too few data is most likely to weaken the sensitivity of the 
search for an energy-dependent anisotropy [|25ll . To fulfill both requirements, the size of the 
energy intervals is chosen to be A log 10 (£') = 0.3 below 8 EeV, so that it is larger than the 
energy resolution even at low energies. At higher energies, to guarantee the determination 
of the amplitude measurement within an uncertainty a - 2%, all events 5,000) with 
energies above 8 EeV are gathered in a single energy interval. 

The amplitude r at the sidereal frequency as a function of the energy is shown in Fig. [51 
together with the corresponding probability P(> r) to get a larger amplitude in each energy 
interval for a statistical fluctuation of isotropy. The dashed line indicates the 99% C.L. 
upper bound on the amplitudes that could result from fluctuations of an isotropic distribu- 
tion. It is apparent that there is no evidence of any significant signal over the whole energy 
range. A global statement refering to the probability with which the 6 observed amplitudes 
could have arisen from an underlying isotropic distribution can be made by comparing the 
measured value K = J^ =1 fc , (where the sum is over all 6 independent energy intervals) 
with that expected from a random distribution for which (K) = 6 [26]. The statistics of 2K 
under the hypothesis of an isotropic sky is ax 2 with 2x6 = 12 degrees of freedom. For 
our data, 2K =19.0 and the associated probability for an equal or larger value arising from 
an isotropic sky is ^ 9%. 

The phase c/> of the first harmonic is shown in Fig. [6] as a function of the energy. While 
the measurements of the amplitudes do not provide any evidence for anisotropy, we note 
that the measurements in adjacent energy intervals suggest a smooth transition between 
a common phase of - 270° in the first two bins below 1 EeV compatible with the 
right ascension of the Galactic Center ace — 268.4°, and another phase (a - 100°) above 



16 



5 EeV. This is intriguing, as the phases are expected to be randomly distributed in case 
of independent samples whose parent distribution is isotropic. Knowing the p.d.f. of phase 
measurements drawn from an isotropic distribution, po((p) = (2n)~ l , and drawn from a 
population of directions having a non-zero amplitude r with a phase (p , pi((p; r , (fo) 
the likelihood functions of any of the hypotheses may be built as : 

i=l i=l 

Without any knowledge of the expected amplitudes r (E) in each bin, the values consid- 
ered in L\ are the measurements performed in each energy interval. For the expected phases 
(fo(E) as a function of energy, we use an arctangent function adjusted on the data as illus- 
trated by the dashed line in Fig. |6] Since the smooth evolution of the phase distribution is 
potentially interesting but observed a posteriori, we aim at testing the fraction of random 
samples whose behaviour in adjacent energy bins would show such a potential interest but 
with no reference to the specific values observed in the data. To do so, we use the method 
of the likelihood ratio test, computing the -2 ln(A) statistic where A = L jh\ . Using only 
Nbins = 6, the asymptotic behaviour of the -2 ln(A) statistic is not reached. Hence, the p.d.f. 
of -2 ln(/l) under the hypothesis of isotropy is built by repeating exactly the same proce- 
dure on a large number of isotropic samples : in each sample, the arctangent parameters 
are left to be optimised, and the corresponding value of -2 ln(A) is calculated. In that way, 
any alignments, smooth evolutions or abrupt transitions of phases in random samples are 
captured and contribute to high values of the -2 ln(A) distribution. The probability that the 
hypothesis of isotropy better reproduces our phase measurements compared to the alter- 
native hypothesis is then calculated by integrating the normalised distribution of -2 ln(/l) 
above the value measured in the data. It is found to be 2 • 10 3 . 

It is important to stress that no confidence level can be built from this report as we 
did not perform an a priori search for a smooth transition in the phase measurements. To 
confirm the detection of a real transition using only the measurements of the phases with an 
independent data set, we need to collect =^1.8 times the number of events analysed here to 
reach an efficiency of - 90% to detect the transition at 99% C.L. (in case the observed effect 
is genuine). It is also worth noting that with a real underlying anisotropy, a consistency of 
the phase measurements in ordered energy intervals is expected with lower statistics than 



the detection of amplitudes significantly standing out of the background noise [|26L 12711 . This 
behaviour was pointed out by Linsley, quoted in [26] : "if the number of events available 
in an experiment is such that the RMS value of r is equal to the true amplitude, then in a 
sequence of experiments r will be significant (say P(> r) < 1%) in one experiment out of 
ten whereas the phase will be within 50° of the true phase in two experiments out of three." 
We have checked this result using Monte Carlo simulations. 

An apparent constancy of phase, even when the significances of the amplitudes are rel- 
atively small, has been noted previously in surveys of measurements made in the range 



10 14 < E < 10 17 eV [28, 29fl. In Ii28ll Greisen and his colleagues comment that most exper- 
iments have been conducted at northern latitudes and therefore the reality of the sidereal 
waves is not yet established. The present measurement is made with events coming largely 
from the southern hemisphere. 
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Figure 7: Rayleigh analysis above 1 EeV using a reduced data set built to get at a constant exposure 
in right-ascension, using on-times of 98% (triangles) and 99% (squares) (see text). Left : analysis 
of the amplitude. The full (dashed) line indicates the 99% C.L. upper bound on the amplitudes 
that could result from fluctuations of an isotropic distribution when using on-times of 98% (99%). 
Right : analysis of the phase, the dashed line being the same as the one plotted in Fig. [6] Results of 
Fig. [6] are also shown with circles. 



5.4. Additional cross-checks against systematic effects above 1 EeV 

It is important to verify that the phase effect is not a manifestation of systematic effects, 
the amplitudes of which are at the level of the background noise. We provide hereafter 
additional studies above 1 EeV, where a few tests can cross-check results presented in 
Fig.® 

The first cross-check is provided by applying the Rayleigh analysis on a reduced data 
set built in such a way that its corresponding exposure in right ascension is uniform. This 
can be achieved by selecting for each sidereal day only events triggering an unitary cell 
whose on-time was almost 100% over the whole sidereal day. To keep a reasonably large 
data set, we present here the results obtained for on-times of 98% and 99%. This allows 
us to use respectively ^77% and 63% of the cumulative data set without applying any 
correction to account for a non-constant exposure. The results are shown in Fig. [7] when 
considering on-time of 98% (triangles) and 99% (squares). Even if more noisy due to the 
reduction of the statistics with respect to the Rayleigh analysis applied on the cumulative 
data set, they are consistent with the weighted Rayleigh analysis and support that results 
presented in Fig. |6] are not dominated by any residual systematics induced by the non- 
uniform exposure. 

From the Fourier analysis presented in section 15.21 we have stressed the decoupling 
between the solar frequency and both the sidereal and anti-sidereal ones thanks to the fre- 
quency resolution reached after 6 years of data taking. However, as the amplitude of an 
eventual sideband effect is proportional to the solar amplitude 112 ill , it remains important 
to estimate the impact of an eventual sideband effect persisting even after the energy cor- 
rections. To probe the magnitude of this sideband effect, we use 10,000 mock data sets 
generated from the real data set (with energies corrected for weather effects) by randomis- 
ing the arrival times but meanwhile keeping both the zenith and the azimuth angles of each 
original event. This procedure guarantees the production of isotropic samples drawn from 
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Figure 8: Top :Anti-sidereal analysis of mock samples built such that a spurious anti-sidereal am- 
plitude stands out from the noise through the sideband mechanism between 1 and 2 EeV (see text). 
Thin histograms : standard Rayleigh analysis of the amplitudes (left panel) and phases (right panel) 
for isotropic samples. Thick histograms : same for samples biased by the energy variations in- 
duced by the atmospheric changes. Dashed histograms : same by reversing the energy corrections. 
Filled histograms : same by amplifying by 10 the energy variations. Bottom : Rayleigh analysis 
of the anti-sidereal amplitude (left) and phase (right) of the real data sample without applying the 
energy corrections (squares), by applying those corrections (triangles), and by reversing the energy 
corrections (circles). 



a uniform exposure with the same detection efficiency conditions than the real data. The 
results of the Rayleigh analysis applied to each mock sample between 1 and 2 EeV at the 
anti-sidereal frequency are shown by the thin histograms in top panels of Fig. [81 displaying 
Rayleigh distributions for the amplitude measurements and uniform distributions for the 
phase measurements. Then, after introducing in to each sample the temporal variations of 
the energies induced by the atmospheric changes according to Eqn. SI it can be seen on the 
same graph (thick histograms) that the amplitude measurements are almost undistinguish- 
able with respect to the reference ones, while the phase measurements start to show to a 
small extent a preferential direction. The same conclusions hold when reversing the energy 
corrections (dashed histograms), but resulting in a phase shift of 180°. Finally, the filled 
histograms are obtained by amplifying by 10 the energy variations induced by the atmo- 
spheric changes. In this latter case, the large increase of the solar amplitude induces a clear 
signal at the anti-sidereal frequency through the sideband mechanism, as evidenced by the 



19 



180 








1 


90 : 








1 


: 


my 

A 


i 




270: 




> 






180 : 








■ Rayleigh analysis 
-•- East/West analysis 



2 3 4 567 10 20 30 



E [EeV] 



300- 



100 



— Anti-Sidereal Frequency 
Solar Frequency 



N=200,000 



200 : -T|; Li ! = : :: 



East/West analysis - 1 0x reversed corrections 

50 100 150 200 250 300 350 

Phase H 



Figure 9: Left : East/West analysis at the sidereal frequency above 1 EeV. Right : East/West anal- 
yses at both the solar (dashed) and the anti-sidereal (thick) frequencies of mock samples built by 
amplifying by 10 the energy variations induced by weather effects. 

distributions of both the amplitudes and the phases. The sharp maximum of the phase dis- 
tribution points towards the spurious direction, while the amplitude distribution follows a 
non-centered Rayleigh distribution with parameter ^ 1.4 x 10" 2 . The spurious mean am- 
plitude stands out from the noise (~ V^/200, 000 - 4 x 10~ 3 ) sufficiently to allow us to 
estimate empirically the original effect to be ten times smaller, at the level of - 1.4 x 10~ 3 . 
This will impact the analyses only in a marginal way even if we had not performed the 
energy corrections, or if we had over-corrected the energies. The energy corrections nec- 
essarily reduce even more the size of the sideband amplitude, well below 10" 3 . Hence, 
only small changes are expected in the anti-sidereal phase measurements on the real data 
when applying (or not) or reversing the energy corrections. This is found to be the case, as 
illustrated in the bottom panels of Fig. [8] In addition, this phase distribution (bottom right 
panel of Fig. [8]) does not show any particular structure aligned with the spurious direction. 
This second cross-check support the hypothesis that the phase measurements presented in 
Fig.[6]are not dominated by any residual systematics induced by the sideband mechanism. 

In the left panel of Fig. [9] we show the results of the East/West analysis (circles). They 
provide further support of the previous analyses (squares) . As previously explained, this 
method relies on the high symmetry between the Eastern and Western sectors. However, 
the array being slightly tilted (^ 0.2°), this symmetry is slightly broken, resulting in a small 
shift between the Eastward and Westward counting rates. As this shift is independent of 
time, it does not impact itself in the estimate of the first harmonic. However, it is worth 
examining the effect of the combination of the tilted array together with the spurious mod- 
ulations induced by weather effects, as this combination may mimic a real East/West first 
harmonic modulation at the solar frequency. The size of such an effect can be probed by 
analysing the mock data sets built by amplifying by 10 the energy variations induced by 
the atmospheric changes. The results obtained between 1 and 2 EeV at both the solar and 
the anti-sidereal frequencies are shown in the right panel of Fig. [9]: while the phase distri- 
bution starts to show a preferential direction at the solar frequency (dashed histogram), the 
same distribution is still uniform at the anti-sidereal one (thick histogram). Hence, it is safe 
to conclude that the results obtained at the sidereal frequency by means of the East/West 
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Table 2: Results of first harmonic analyses in different energy intervals, using the East/West analysis 
below 1 EeV and the Rayleigh analysis above 1 EeV. 



method are not affected by any systematics. 

5.5. Results at the sidereal frequency in cumulative energy bins 

Performing the same analysis in terms of energy thresholds may be convenient for opti- 
mizing the detection of an eventual genuine signal spread over a large energy range, avoid- 
ing the arbitrary choice of a bin size Alog 10 (£'). The bins are however strongly correlated, 
preventing a straightforward interpretation of the evolution of the points with energy. The 
results on the amplitudes are shown in Fig. \\0\ They do not provide any further evidence 
in favor of a significant amplitude. 

6. Upper limits and discussion 

From the analyses reported in the previous Section, upper limits on amplitudes at 99% 
C.L. can be derived according to the distribution drawn from a population characterised by 
an anisotropy of unknown amplitude and phase as derived by Linsley [fl3ll : 

[2 1 r ruL As T (rs\ ( s 2 + r 2 /2\ „ r 

where 7 is the modified Bessel function of the first kind with order 0, and cr = ^2 IN 
in case of the Rayleigh analysis, and cr EW = (n(cos5) /2(sin#)) x ^2/N in case of the 
East/West analysis. 

As discussed in the Appendix, the Rayleigh amplitude measured by an observatory 
depends on its latitude and on the range of zenith angles considered. The measured am- 
plitude can be related to a real equatorial dipole component d ± by d ± - r/(cosS). This 
is the physical quantity of interest to compare results from different experiments and from 
model predictions. The upper limits on d , are given in Tab. |2] and shown in Fig{TT] together 
with previous results from EAS-TOP O, KASCADE KASCADE-Grande [Q and 
AGASA pih . and with some predictions for the anisotropics arising from models of both 
galactic and extragalactic UHECR origin. The results obtained in this study are not consis- 
tent with the - 4% anisotropy reported by AGASA in the energy range 1 < E/EeV < 2. 

If the galactic/extragalactic transition occurs at the ankle energy QlJ], UHECRs at 1 EeV 
are predominantly of galactic origin and their escape from the galaxy by diffusion and 



21 



CD 
"O 
=i 

"Q. 

E 
< 



Rayleigh Analysis 
East/West Analysis 



10 



1 _ 



-r 



10 



2 _ 



-f 



10 



-3 . 



0.3 



t i 



2 3 4 56 10 20 
E th [EeV] 



O 



10" 1 



10" : 



c 180 

CD 
C/) 
03 

£ 90 



Rayleigh Analysis 
East/West Analysis 



0.3 



1 2 3 4 56 10 20 

E th [EeV] 



Rayleigh Analysis 
East/West Analysis 



270 



180 



T T f 



• T • 



0.3 



1 2 3 4 56 10 20 

E th [EeV] 



Figure 10: Same as Fig.[5]and Fig. [6] but as a function of energy thresholds. 
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Figure 1 1 : Upper limits on the anisotropy amplitude of first harmonic as a function of energy from 
this analysis. Results from EAS-TOP, AGASA, KASCADE and KASCADE-Grande experiments 
are displayed too. An analysis of the KASCADE-Grande data with the East/West method delivers an 
additional limit for 3 10 15 eV. Also shown are the predictions up to 1 EeV from two different galactic 
magnetic field models with different symmetries (A and S), the predictions for a purely galactic 
origin of UHECRs up to a few tens of 10 19 eV (Gal), and the expectations from the Comp ton- 
Getting effect for an extragalactic component isotropic in the CMB rest frame (C-G Xgal). 



drift motions are expected to induce a modulation in this energy range. These predictions 
depend on the assumed galactic magnetic field model as well as on the source distribu- 
tion and the composition of the uhecrH. Two alternative models are displayed in Fig. 
ITTl corresponding to different geometries of the halo magnetic fields H. The bounds re- 
ported here already exclude the particular model with an antisymmetric halo magnetic field 
(A) and are starting to become sensitive to the predictions of the model with a symmetric 
field (5). We note that those models assume a predominantly heavy composition galactic 
component at EeV energies, while scenarios in which galactic protons dominate at those 
energies would typically predict anisotropics larger than the bounds obtained in Fig. ITTl 
Maintaining the amplitudes of such anisotropics within our bounds necessarily translates 
into constraints upon the description of the halo magnetic fields and/or the spatial source 
distribution. This is particularly interesting in the view of our composition measurements 
at those energies compatible with a light composition 113411 . Aternatively to a leaky galaxy 
model, there is still the possibility that a large scale magnetic field retains all particles in 



5 The dependence of the detection efficiency on the primary mass below 3 EeV could affect the details of 
a direct comparison with a model based on a mixed composition. 
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the galaxy 11351 . 13611 . If the structure of the magnetic fields in the halo is such that the turbu- 
lent component predominates over the regular one, purely diffusion motions may confine 
light elements of galactic origin up to ^ 1 EeV and may induce an ankle feature due to the 
longer confinement of heavier elements at higher energies QUO - Typical signatures of such 
a scenario in terms of large scale anisotropics are also shown in Fig. ITTI ("dotted line) : the 
corresponding amplitudes are challenged by our current sensitivity. 

On the other hand, if the transition is taking place at lower energies around the sec- 
ond knee at 5 x 10 17 eV [7], UHECRs above 1 EeV are dominantly of extragalactic 
origin and their large scale distribution could be influenced by the relative motion of the 
observer with respect to the frame of the sources. If the frame in which the UHECR distri- 
bution is isotropic coincides with the CMB rest frame, a small anisotropy is expected due 
to the Compton-Getting effect. Neglecting the effects of the galactic magnetic field, this 
anisotropy would be a dipolar pattern pointing in the direction a - 168° with an amplitude 
of about 0.6% 113811 . On the contrary, when accounting for the galactic magnetic field, this 
dipolar anisotropy is expected to also affect higher order multipoles 113 911 . These ampli- 
tudes are close to the upper limits set in this analysis, and the statistics required to detect 
an amplitude of 0.6% at 99% C.L. is — 3 times the present one. 

Continued scrutiny of the large scale distribution of arrival directions of UHECRs as a 
function of energy with the increased statistics provided by the Pierre Auger Observatory, 
above a few times 10 17 eV, will help to discriminate between a predominantly galactic or 
extragalactic origin of UHECRs as a function of the energy, and so benefit the search for the 
galactic/extragalactic transition. Future work will profit from the lower energy threshold 
that is now available at the Pierre Auger Observatory [40J. 



This article is dedicated to Gianni Navarra, who has been deeply involved in this study 
for many years and who has inspired several of the analyses described in this paper. His 
legacy lives on. 

Appendix 

The first harmonic amplitude of the distribution in right ascension of the detected cos- 
mic rays can be directly related to the amplitude d of a dipolar distribution of the cosmic 
ray flux of the form J (a, 8) = 7o(l + d d ■ u), where u and d denote respectively the unit 
vector in the direction of an arrival direction and in the direction of the dipole. Eqn.QTJand 
Eqn. [12] can be used to express a, b and N as : 

2 r^max nix 

a = — I d5 da cos 8 J(a, 8) to{8) cos a, (19) 
N Js rni „ Jo 

= — / dd J da cos 8 J(a, 8) to{8) sin a, 
N J<s„„„ Jo 



u iain 



to 

d8 J da cos 8 J {a, 8) a>(8) 

Smin JO 
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where we have here neglected in the exposure co the small dependence on right ascension. 
Writing the angular dependence in J{a, 5) as d ■ u = cos 6 cos 8 d cos(a - a d ) + sin 5 sin 5 d , 
with 6 d the dipole declination and a d its right ascension, and performing the integration in 
a in the previous equations, it can be seen that 



Ad ± 



(20) 



1 + Bd z 
where 

f dd a>(6) cos 2 6 f d5 a>(6) cos 5 sin 5 

A = -r > B = —r 

J dd a>(6) cos 6 J dd co{5) cos 6 

and d z = d sin 8 d denotes the component of the dipole along the Earth rotation axis while 
d ± = dcos6 d is the component in the equatorial plane fl30h . The coefficients A and B 
can be estimated from the data as the mean values of the cosine and the sine of the event 
declinations. In our case, A = (cos 6) - 0.78 and B = (sin 6) - -0.45. For a dipole 
amplitude d, the measured amplitude of the first harmonic in right ascension r thus depends 
on the region of the sky observed, which is essentially a function of the latitude of the 
observatory t site , and the range of zenith angles considered. In the case of a small Bd z 
factor, the dipole component in the equatorial plane d ± is obtained as d ± - r/ (cos 6). The 
phase ip corresponds to the right ascension of the dipole direction a d . 

Turning now to the East-West method, the measured flux from the East sector for a 
local sidereal time a can be similarly expressed as 

I E (a°)= d(f> d6sm6e{e)J(e,(p,a\ (21) 
J-tt/2 Jo 

and analogously for the measured flux coming from the west sector changing the azimuthal 
integration to the interval [n/2, 3n/2]. Expressing d ■ u in local coordinates (6, (p and or ), 
and performing the integration over (p we obtain for the leading order 

h - Iw , (k 2d ± C o 

(« ) = — t: -— — — - — - sin((r - a d ), (22) 



(1 E +I W ) MX + d z D sin € site ) 

where 

f d9 e(9) sin 2 9 [ d9 e(6) sin 6 cos 9 

C = - r , D = - — . 

J dO e(9) sin 9 J d9 e{9) sin 9 

In this calculation any dependence of the exposure on the local sidereal time a gives 
at first order the same contribution to the East and West sectors flux, and thus gives a 
negligible contribution to the flux difference^. The next leading order, proportional to the 
equatorial dipole component times the sidereal modulation of the exposure, is negligible. 
The coefficients C and D can be estimated from the observed zenith angles of the events. 



6 The exposure dependence on azimuth present in Auger due to trigger effects at low energies has a 2n/6 
frequency that makes it cancel out in the computation of the first harmonic, and thus does not affect the above 
result. 
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In our case, C = (sin#) 0.58 and D = (cos 6) - 0.78. The total detected flux averaged 
over the local sidereal time can be estimated as (I E + I w ) = N/2n. In case Dd z <sc 1, we get 
finally : 

(I E -I w )(a^ = -P-^^sm(a°-a d ). (23) 
2n n 
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